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We study the electronic structure and the phase diagram of non-interacting fermions confined to 

00 ' hexagonal optical lattices. In the first part, we compare the properties of Dirac points arising in 
, the eigenspectrum of either honeycomb or triangular lattices. Numerical results are complemented 

■ by analytical equations for weak and strong confinements. In the second part we discuss the phase 
CN ■ dia gram and the evolution of Dirac points in honeycomb lattices applying a tight-binding description 

^ I with arbitrary nearest-neighbor hoppings. With increasing asymmetry between the hoppings the 

^ , Dirac points approach each other. At a critical asymmetry the Dirac points merge to open an energy 

gap, thus changing the topology of the eigenspectrum. We analyze the trajectory of the Dirac points 
and study the density of states in the different phases. Manifestations of the phase transition in the 
, temperature dependence of the specific heat and in the structure factor are discussed. 

r-^ ! I. INTRODUCTION 

. . . , . . . . 

^ • The isolation of a single layer of graphene , which is a two-dimensional honeycomb lattice of carbon atoms, has at- 

' tracted considerable attention, both for its basic interest and because it may pave the way to carbon-based electronics. 
Cj ■ The low energy electronic properties are described by the massless Dirac equation, which causes many anomalies with 
. ' respect to semiconductor physics^"^. Close to half-filling, the band structure near the Fermi energy is given by two 
"j^ , degenerate Dirac cones, and the centers of these cones (called Dirac points) are located at the two distinct K-points 
' of the first Brillouin zone. A first experimental validation of the relativistic band structure was the measurement of 
I I the anomalous sequencing of the Quantum Hall plateaus^'^. 
'■^ ' In graphene there is little room to tune system parameters such as the interaction strength or the hopping am- 
^ [ plitudes, which prevents a systematic study of its phase diagram. In contrast for cold fermionic atoms confined in 
O , honeycomb optical lattices'" the parameters are highly tunable and many of the theoretically predicted phases might 

1 be realized there. Recent theoretical works on honeycomb lattices have studied a topological phase transition in 
the single particle spectrum which is due either to asymmetric hopping energies®, or to a Kekule distortion of the 
hoppings^'^", or to the application of an ac electric field^^. Further work on honeycomb optical lattices has dealt with 
the realization of the anomalous Hall and the Spin Hall effects^^, the Wigner crystallization of Tr-bands^^'^**, or with 

■ the Quantum Hall effect without Landau levels^^"^^. Finally, we note recent theoretical work on the appearance and 



' manifestation of Dirac cones in the band structure of triangular photonic crystals^®^^^. 

There are two main purposes of this work. The first one is to compare the Dirac points present in the lower lying 
bands of triangular and honeycomb lattices. Therefore we first discuss in section II the potential landscape of these 
lattices. Then, in section HI we discuss the Dirac cones appearing in both lattice types and derive analytical results 
within the nearly-free-particle and tight-binding limits. The second main goal of this work is to analyze in detail the 
, topological phase transition in honeycomb lattices with asymmetric hopping energies®. In section IV, the evolution 
of the Dirac points and the resulting phase diagram for the honeycomb lattice as a function of the relative tunnel 
couplings are studied within the tight-binding limit . The density of states for the different phases is derived in section 
V. Experimental signatures of the phase transition such as the structure factor and the temperature dependence of 
5— i ' the specific heat are explained in section VI. Finally, we discuss the topological structure of the phase transition in 
section VII. Some technical derivations are left to the Appendices. 
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II. POTENTIAL LANDSCAPES 



We consider the optical lattice that was experimentally realized by Grynberg et al7 . It is created by three interfering 
traveling laser beams. The corresponding wave vectors are in the x — y plane and form equal angles between them, 
as illustrated in Fig. la. For simplicity we assume the polarization of the beams to be perpendicular to the plane, 
although the confining potential is independent of the polarization for large enough detuning^^. 
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FIG. 1: (a) Three traveling lasers with wavevectors interfere and build up the hexagonal lattice with reciprocal unit vectors 
bi,b2, and bs = — (bi + b2). (b) Unit cell in real space spanned by the unit vectors ai,a2. For triangular lattices the minima 
are located at the center of the hexagons, while for symmetric honeycomb lattices minima lie at the A and B sites, (c) Unit 
cell in reciprocal space (first Brillouin zone) indicating the position of the symmetry points F, Mi, M2, M3 and the K-points. 



The amplitude of the total electric field (polarized in the z-dircction) and the light intensity can then be written 

E{r) = ^i?,exp(zq, -r) (1) 
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I{r) = \E{r)\^ ^ El + El + eI + 2E1E2 cos(b3 • r) + 2E2E:i cos(bi • r) + 2E3E1 cos(b2 • r) . (2) 

As illustrated in Fig. la bi = q2 — qa plus cyclic permutations. We note that the system is essentially insensitive to 
phase drifts between the lasers, since these just shift the lattice without changing the potential landscape. We can 
therefore assume Ei> Q without loss of generality. 

The optical confining potential ^(r) is proportional to the intensity of the laser field and the proportionality factor 
is positive (negative) if the laser frequencies are blue (red) detuned with respect to the transition frequency of the 
atoms. Up to a constant the confinement is given by 

V{v) = ^F,cos(b,-r), (3) 


with Vi oc -E2-E'3 plus cyclic permutations and > (V,- < 0) for blue (red) detuning. 

The potential V{t) has an underlying triangular Bravais lattice, whose lattice spacing a is determined by the 
wavelength of the traveling lasers a — An/Sq = 2A/3. The unit vectors for real and reciprocal lattices are given by: 

-K;^)^ *"^7s(f)^ "-M-f) 

with bi • SLj = 27r(5y . The unit cells of both reciprocal and real space can be chosen of honeycomb form as shown in 
Figs lb and c. 

In the case of red detuning, the potential minima lie at the centers of the hexagons thus forming a triangular 
lattice for all values of Vj < 0. If the lasers are blue detuned and their intensities are equal, then the minima of 
the potential V{r) are located at the vertices of the honeycombs while there is a maximum of the potential at the 
center of each honeycomb, see Fig. lb. If the intensities Vj become different (i.e. the electric field amplitudes Ej 
begin to be different), the two minima per unit cell approach each other until they merge for strong asymmetries, 
namely, when \Ei — i?2| > E3 01 E3 > Ei + E2- Thus the potential V{r) (3) has two minima per unit cell for 
\Ei ~ E2\ < Es < El + E2- It is evident from Eq. (3) that the potential ^(r) has inversion symmetry, which is crucial 
for the stability of the Dirac points^^. The position of the minima are given by r = irmin with: 



/ ^ sgn(£;2 - El) arccos ( [mEMEM^\ \ 



E3 V 4B1B2 



— arccos I -y 



(5) 



We will show in section IV that the energy dispersion i?k calculated within the tight-binding approximation has the 
same structure as the potential V{y). The analog of the motion of the potential minima for the case of the energy 
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(a) V = 4 (b) V = 0.4 (c) V = 0. (d) V = -0.4 
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FIG. 2: Band structure for the symmetric lattice (Vi = V2 = V3 = V) along characteristic lines of the first Brillouin zone for 
various values of V. Symmetry points are explained in Fig. Ic. The energy of Dirac point is set to and the unit of energy is 
(ftV2m)(47rVa^). 



spectrum is the evolution of the Dirac points (8), which will finally merge, causing a topological phase transition from 
a semimetal to an insulator. 



III. SYMMETRIC LATTICES 



Because of periodic translational invariance, the crystal momentum is conserved and the eigenspectrum of a particle 
moving in the potential V{r) of Eq. (3) is most conveniently calculated in Fourier space. The eigenfunctions are thus 
written as V^qCr) = J2m n '^q-Gmn 6xp(i(q — Gmn) ■ r), where q is restricted to be in the first Brillouin zone, Gmn = 
mbi + nh2 denotes a general reciprocal lattice vector, and the potential is written as V{y) = X^mn ^mn exp(iG„i„ • r) 
where = i [Vz{5m,i5n,i + 5m-i5n-i) + Vi5nfi{5m,i + 5m -1) + V25m,o{5n,i + 5n-i)] are the Fourier components 
of the potential. Due to inversion symmetry of the lattice Vmn = V-m.-n- 

For each q the Fourier components Cq-G„„ form a vector Vq, where each entry is labeled by a specific pair of 
integers (vq)m„ = Cq_G„„- The Schrodinger equation then reduces to a set of linear equations for each Vq given by: 

MqVq = i;qVq (6) 

(-'^q)(mi,ni),(m2,n2) ~ '^mi ,"12 '^"l ,"2 -^q— Gm ^ _„ ^ + ^^i— m2,"l— "2 

where = h?q^/2m denotes the kinetic energy. 

Figure 2 shows the band structure for symmetric potentials Vi = V2 = V3 = V for various confinement strengths 
V. Most importantly, we note that, both for triangular and honeycomb lattices, some pairs of bands touch at the 
K-points and in the vicinity of those touching points the band structure is given by cones, which are called Dirac 
cones because of the similarity to the relativistic energy dispersion. 

There are however important differences between both cases. For the honeycomb lattice the first pair of Dirac 
points (located at ±K) arises due to the touching of the two lowest bands. Furthermore, the Fermi surface in the 
relevant energy interval is exclusively determined by these Dirac cones, leading in particular to a vanishing density of 
states at complete filling of the lowest band. 

By contrast, for a triangular lattice the first pair of Dirac points is caused by the touching between the second 
and third band. Since these bands are also degenerate at the F point, the Dirac points are not isolated but resonate 
with a continuous band that dominates the density of states in the energy regime of the Dirac cones. Furthermore, 
the filling factor (number of atoms per unit cell) at which the Fermi energy coincides with the Dirac point is not an 
integer, also in contrast to the honeycomb lattice. 

The formation of the Dirac cone can be derived within the nearly-free-particle approximation, where the potential 
is treated perturbatively as shown in Appendix A. The distinction between honeycomb and triangular lattices is 
however better appreciated in the tight-binding regime, where the Hamiltonian is written in terms of Wannier states, 
each localized at a potential minimum as shown in Appendix B. 

In the following we concentrate on the honeycomb lattice within the tight-binding approximation. 
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FIG. 3: (a) Phase diagram of the asymmetric honeycomb lattice within tight-binding description. Black lines with arrows are 
coordinate axes for the tunneling amplitudes. Pink central region denotes semimetal with two distinct Dirac points. White 
regions labeled by Ai,A2, A3 denote insulating phases, in which no Dirac points are present. In the insulating phase there is a 
single band minimum in the first Brillouin zone located at Mi for Ai and so forth, (b) Motion of Dirac points as a function of 
ts for ti/t2 = 1 (red line in a and b) ti/t2 =3/2 (blue) ti/t2 ~ 2/3 (green). Arrows in (a) and (b) indicate direction of motion 
of Dirac points for increasing t^. The symmetric case ti — t2 — ts (like in graphene) is indicated by red circles. 



IV. DIRAC POINTS IN ASYMMETRIC HONEYCOMB LATTICES 



We now discuss the eigenspectrum of a honeycomb lattice with asymmetric hoppings, described by the following 
tight-binding Hamiltonian: 



3 

where Ok, ^k destroy Bloch waves on the two different sublattices and tj > denote nearest-neighbor hopping energies. 
For symmetric hoppings the standard Hamiltonian (Bl) is rediscovered. 

A possible realization of this Hamiltonian is a deep optical potential ^(r) (3) with different values of Vj. Due to 
the inversion symmetry of the confinement potential (3) the onsite energies at the two sublattices are equal and thus 
can be set to zero in Eq. (7). For asymmetric laser intensities the potentials barriers separating a given potential 
minimum from the three nearest minima are different. In particular, the corresponding distances between neighboring 
minima are also different, as follows from Eq. (5). Both effects will lead to asymmetric nearest-neighbor hoppings. 
We note that a pure shift in the position of the potential minima without a change in the magnitude of the tunnel 
hopping energies only gives rise to a phase factor in ^k (7) and thus does not affect the energy spectrum £^k = i|'?fc|- 

The energy spectrum of the Hamiltonian (7) is symmetric around zero energy and is given by = ilCfcl- We note 
the close analogy between the energy amplitude ^k and the total electric field E{r). This causes a close similarity 
between = |^kP and V{r) oc E{r)'^. As discussed above for the symmetric case ti = t2 = ^3, valence and conduction 
band (defined by negative and positive energies, respectively) touch at the two K-points, where they form Dirac cones. 
Introducing an asymmetry between the tunneling amplitudes the two Dirac points move away from the K-points and 
approach each other. At a critical asymmetry the Dirac points merge at one of the three inequivalent M-points, and 
by increasing the asymmetry further a band gap opens. 

Figure 3 illustrates the phase diagram and the displacement of the Dirac points as a function of the tunnel ampli- 
tudes. The central pink region in Fig. 3a embraces the parameter regime of the semimetallic phase, where the asymme- 
try is small enough for the two inequivalent Dirac points to exist. It is defined by the condition |ti — < < ti+t2, 
which remains invariant under the permutation of the three subindices. The Dirac points are located at k = ik^i, 
with: 



/ 2 f /t?-(t2-tl)A \ 

I - arccos I — w ^ ^ — — I 



kD = 



4tit 



11:2 



^ ;^sgn(^,-^2) arccos (i^v/^^^ 



(8) 
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FIG. 4: Density of states per unit cell for different hoppings normalized to ti + + = 1- (a) Linear energy dependence at 
low energy within semimetallic phase, (b) Square root energy dependence at phase transition ts = 0.5. (c) Insulating phase. 
Blue curves are analytical approximations to red curves given in main text. 



If necessary the wavevector k£) as defined above is supposed to be folded back to the first Brillouin zone by adding a 
uniquely defined reciprocal lattice vector. 

Figure 3b shows the displacement of the Dirac points following three different lines of the phase diagram, each 
defined by a fixed ratio ti/i2 while spanned by varying t^. We focus first on the red line, which corresponds to the 
partially symmetric case ti = t2 previously studied in Ref. [8]. For ^3 = valence and conduction band touch along the 
extended dashed line. Increasing ^3 causes the two inequivalent Dirac cones to approach each other along horizontal 
lines and for symmetric hoppings ^1=^2= ^3 (center of triangle in Fig. 3a) they are located at the K-points. Finally 
the Dirac cones merge at the symmetry point M3 = (bi + b2)/2 when ^3 = ^1+^2- 

In the generic case of ii 7^ ^2 7^ ^3 one finds the motion shown by the blue and green lines. We now discuss the 
motion for increasing t^^. For ^3 < \ti — ^2! the system is an insulator and there is a single band minimum located at 
Ml, if ti > t2 (blue line in region Ai of Fig. 3a), or M2, if ^2 > ti (green line in region A2 of the same figure). While 
in regions Ai or A2, the minimum stays pinned at the fixed points Mi or M2 and motion along the blue or green 
line only causes the corresponding gap to decrease until it becomes zero when the central pink triangle is touched 
(phase transition). For \ti — ^2! < ^3 < +^2, in both cases (blue and green) the system enters the central triangle of 
Fig. 3a defining the semimetallic phase, where the band minimum at Mi or M2 separates into two Dirac points. The 
position of the Dirac points are related to each other by the inversion symmetry around the symmetry points F and 
Ml, M2, M3. Finally at ^3 =^1+^2, where the lines leave the pink triangle in order to enter region A3, the Dirac 
points merge at M3 and the system undergoes a phase transition from a semimetal to an insulator. Further increase 
of ^3 only causes the gap at M3 to increase. 



V. DENSITY OF STATES 



We now discuss the density of states, which is defined by p{E) — A " -^k), where A = NAu is the 

area, N the number of unit cells and A^ = a^-\/3/2 the size of the unit cell. Characteristic energies are given by the 
eigenenergies at the symmetry points i^Ma — 1^3 — ii— ^2!, ^-Mi = l^i — ^2 — i3|,£'M2 = 1^2 — ^1 — ^31 and i^r = ^1+^2 + ^3- 
Ey is the energy maximum and thus sets the bandwidth. The density has a step at _E = ztE'r- 

Without loss of generality we focus on the case > ti,t2, so that the energy always has a saddle point at the 
symmetry points M2 , M3 , which results in logarithmic van Hove singularities at j Ems ■ Whether M3 is a saddle 
point or a minimum depends on the tunnel amplitudes. For < ti + t2 the energy has a saddle point at M3 and 
there are two Dirac points in the first Brillouin zone. For symmetric hoppings ^3 = ^1= ^2 the band structure close 
to zero energy is given by rotationally symmetric cones centered at the Dirac points. As the hopping energies become 
different from each other, those cones are stretched and the Dirac points move, but still the density of states vanishes 
linearly close to zero energy {E <C Eyi^, Eyi2)- In particular for ti = t2 > H/2 the low energy spectrum is given 

by i?kD+q — ij^^^lq^'+v^q^, Vx — a^/ At\ — i|/2, Vy = at3\/3/2, resulting in a density of states p{E) — gyE/2nVxVy, 

where the degeneracy factor gy = 2 takes into account the two Dirac points. We note that the velocities Vx,Vy are 
different if the hoppings are different. While Vy increases with increasing t^, decreases and cancels at the phase 
transition = 2ti — 2^2- Examples of the density of states in the semimetallic phase are given in Fig. 4a. 

In the opposite limit > ti + t2 the energy has its minimum at M3 and the system is gapped by 2A where 
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A = Ei 



M3 



^3 — ti — 12- Correspondingly the density of states is zero for E < A and has a finite step at = A. For 



ti = t2 <tz/2 the low energy spectrum is given by iJko+q = (A^ 



-vlql+vlqlf/^, = a^JAti/2, Vy = 
leading to a density of states p{E) = 9{E — A)E /2'KVxVy, where 6{x) denotes the Heaviside step function. This 
situation is depicted in Fig. 4c. 

Right at ty, =ti+ 12 the two Dirac cones merge at M3 with Ems = 0, so that M3 is the minimum of the conduction 
band, but there is still no band gap. In the particular case ti=t2= H/2, the density of states can be mostly derived 
analytically: 



^1 



Mg+q — 



p{E) 



2 



6 + 2 cos{qxa) — 8 cos{qxa/2) cos(\/3aq'y/2) 



1/2 



4Vi 



(9) 
(10) 



E I2t-j, denoting the energy in units of the bandwidth. This situation is depicted in Fig. 4b. We note that 
t t-i\i>{aqyl2Y + (a/2)^(q'^ — Gg^g^ — 3g^)/4]^/2, which shows that the confinement in the x-direction (along 
which the Dirac points have merged) is much weaker than in the y-direction. Equation (10) is valid for the whole 
energy range < e < 1. An expansion around = results in: 



with e = 



4v^ 
A„t37r2 



dx 



- 4V^r(5/4) 

^f3A„^2r(3/4) 




(11) 



Since the square root has a diverging slope at the origin, it can be viewed as a transient behavior between the linear 
density of states of the semimetallic phase and the step found for the insulating phase. We note that for e = 0.5, 



which is the energy of the two saddle points _Emi = E\[,^ , the integrand of Eq. (11) contains the factor (1 - 
results in a logarithmic divergence in the density of states. 
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VI. SPECIFIC HEAT AND STRUCTURE FACTOR 



Finally we discuss two experimental signatures of the phase transition expected when fermions populate the lattice 
at half filling, namely, the structure factor accessible by Bragg spectroscopy* and the temperature dependence of the 
specific heat^^. The specific heat is defined by cy = du/dT = J dEp{E)Ef{E), where f{E) denotes the Fermi 
distribution, u the energy density and p{E) the density of states as calculated in the previous section. We limit 
ourselves to the case of half-filling, where the chemical potential p = fov all temperatures because of the symmetry 
of the energy around E = inherent to the employed tight-binding approximation. 

The resulting equations are: 



/ 



dEp{E)Eta.nh{E/kBT) ; cy 



1 



7° 

Jo 



dE- 



p{E)E^ 



(12) 



where we have used p{E) 
difl[erent phases: 



2kBT'^ Jo "'^ cosh^{E/2kBT) 
p{—E). This results in a different temperature dependence of the specific heat for the 



Cv 



T'^ 9gyk%C{3)/2iTVxVy 

{kBT/hf^ 2.1kBlAu 

^ exp{-A/kBT) (6T3 + GT^A + STA^ + A^) 2k%/2-KVxVyT 



for p{E) = gyE/2-KVxVy 
for p(i;)~(E/2i3)V2o.53/A„t3 
for p{E) ~ e{E - A)E/2-KVxVy 



(13) 



where C, denotes the Riemann Zeta Function. Interaction induced corrections to the specific heat in the semimetallic 
phase have recently been studied^^. 

As pointed out in Ref. 8 another way to experimentally visualize the phase transition is to measure the structure 
factor by Bragg spectroscopy. Wc note that the structure factor is up to a constant given by the imaginary part 
of the dynamical susceptibility which is well known for graphene^'^^^e ^^^^^ within the semimetallic phase, the 
results for the susceptibility of graphene can be easily adapted to include asymmetric hoppings. Noting that E^^ ~ 



(u: 



.1/2 



'4ly) ' = "^^^^ = qxVxjKy = qyVy,K 

susceptibility of arbitrary doping by setting vp = ^ and replacing q 
As an illustration, we state the result for half-filling 



KyY^^, one can reuse all formulae of the dynamic 



5(q,w) 



2 ? 
V a 



y^y 



167r 



VxVy. (Jj 



7 9 



9 2 
"yly 



2 2 
VxQx 



Vyqy) 



(14) 
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Here g = 2 is due to the two Dirac points (valley degeneracy). For ti = t2 > H/'i, the velocities are given by 
Vx = a (Ati — tf)^^^ /2, Vy = at3\/3/2. Wc note that in Ref. 8 only the energy spectrum was considered while the 
overlap factor (k, q) resulting from the spinor eigenfunctions was ignored. Therefore the term present in the 
numerator of Eq. (8) of Ref. 8 is absent in Eq. (14) above. We note that ^(k, q) = Iw^+q^A'^kAl^ describes how 
easily the spinor wavcfimction WkA (A = =Fl labels valence and conduction band) can be scattered in to Uk+q.A' and 
is responsible for example for the absence of backscattering and the Klein paradox in graphene.^'' The general form 
of the spinor wavefunction for the tight-binding Hamiltonian (7), is WkA = 2~^/^(— A^k/|^k|, !)• 

For the gaped case ti = t2 < ^3/2, an analytical formula for the structure factor at half-filling and for energies close 

to the gap ^ (A^ + vlql + vlql) ~ A + ^ql + ^ql is given by: 

5(q,a;) = ^(w - 2A - ^ - ^)— — (15) 

While the step function originates from the single particle cigcnspcctrum, the whole q-dependence arises from the 
spinor character of the wavefunctions captured in the factor ^(k, q) described above. 



VII. TOPOLOGICAL PROPERTIES 



Interestingly, a gap only opens after the two Dirac points have merged. Mathematically it is straightforward to 
show that an energy minimum of the conduction band that is not located at one of the symmetry points Mi, M2, 
M3, r is automatically a Dirac point. To show this, wc note that a minimum of the conduction band is a root of 
Vk-E^ = Vk (»?k?/k)' ■^itli ??k = is + ti exp(ikai) + 12 exp(«ka2). 

Vk (??kr?k) = ViRe(?7k) + V2lm(?7k) (16) 
Vl = 2VkRc(r7k) = -2 [ti sin(ai • k)ai + t2 sin(a2 • k)a2] 
V2 = 2VkIm(?7k) = 2 [ti cos(ai • k)ai + t2 cos(a2 • k)a2] 

Since ai,a2 are linear independent, vi,V2 are also linear independent, unless sin((a2 — ai) • k)) = 0, a condition is 
only fulfilled at the symmetry points Mi,M2,M3,r. This analysis shows that a minimum which is not located at 
a symmetry point, is necessarily a Dirac point, since then both the imaginary and the real parts of 77k vanish and 
therefore the energies are i?ip,k = T|??k| =0. Furthermore, these minima always occur in pairs due to time reversal 
symmetry (inversion symmetry around the F-point). 

The deeper reason for the stability of the Dirac points is the momentum space topology of the Bloch 
wavefunctions^"'^'^^. This is nicely illustrated by noting that the Berry phase, defined as^^ (t)B = i (ik(uk| Vkj^k) 
with S denoting a closed path in reciprocal space, takes values ±7r, depending on whether S encloses one Dirac 
point, or the other, or both. Thus each Dirac point is a source of a ±7r delta-function flux of Berry curvature, and 
this flux remains invariant under perturbations that preserve space and time inversion symmetry^^'^^. Only if the two 
Dirac points merge the Berry curvature vanishes within the whole first Brillouin zone and a gap can open. 

We note that the Dirac points are unstable against perturbations that break spatial invariance. Examples are 
substrate-induced gap opening in graphene due to breaking of sublattice symmetry^° or honeycomb lattices where 
nearest-neighbor hoppings are periodically modified to form a Kekule pattern^. 



VIII. SUMMARY AND CONCLUSIONS 



We have studied the electronic structure and the phase diagram of non-interacting fermions confined to hexagonal 
optical lattices. In the first part of the paper, we have analyzed the appearance of Dirac points in the band structure of 
fermionic atoms populating triangular and honeycomb optical lattices for different confinement strengths. Numerical 
results were complemented by analytical equations both for strong and weak potentials. The Dirac points arising in 
honeycomb and triangular lattices differ in several important aspects. While in honeycomb lattice the Dirac cone is 
isolated, it resonates with a continuum band for a triangular lattice. Furthermore, in honeycomb lattices the Fermi- 
energy coincides with the Dirac point for integer filling, while a non-integer filling is needed in the case of a triangular 
lattice. For the honeycomb lattice the first pair of Dirac points arises from the touching of the two lowest bands, 
which within a tight-binding description stem from the ground state of the two potential minima per unit cell. By 
contrast, for a triangular lattice the first pair of Dirac points occurs at the touching of the second and third band, 
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which within a tight-binding description are formed by the doubly degenerate first excited state of the single potential 
minimum per unit cell. 

In the second part of this work, we have focused on the phase diagram and the evolution of the Dirac points in 
honeycomb lattices adopting a tight-binding description with asymmetric nearest-neighbor hoppings. The semimetallic 
phase is not only realized for strictly symmetric hoppings but rather for an extended region of the hopping parameter 
space. This stability of the semimetallic phase is based on topology. With increasing asymmetry the Dirac points 
approach each other along continuous trajectories and at a critical asymmetry they finally merge and an energy gap 
opens. We derive analytic formulae for the trajectories of the Dirac points as well as for the density of states. Right at 
the phase transition the density of states increases like a/B, which interpolates between the linear dependence of the 
semimetallic phase and the step-like increase of the insulating phase. We also show how the phase transition becomes 
manifest in a change of the temperature dependence of the specific heat as well as in the properties of the structure 
factor. 

An interesting follow-up of this work concerns the influence of interactions on the phase diagram, since optical 
lattices permit a systematic control of the effective interaction strength by using Feshbach resonances or modifying 
the lattice properties'^. 
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APPENDIX A: DIRAC CONES IN WEAK SYMMETRIC LATTICES 



Here we apply the nearly-free-particle approximation to derive the eigenspectrum determined by Eq. (6) with 

Vi = V2 = V3 = V close to the K-point. In absence of any potential the dispersion relation at the K-point is three- fold 
degenerate, where the three states can be labeled by the k-values k e {K, K — do, K — Gn}. Treating the lattice 
within lowest order perturbation theory, the eigenvalue problem at k = K -|- q has then the following form: 

E^+^ V/2 V/2 \ ( 

y/2 i^&-G,o+q v/2 CK-Gio+q = i^K+q CK-G.o+q (Al) 

V/2 V/2 VcK-Gii+q/ VcK-G.i+q/ 

The solutions of the above eigenproblem for q = are given by = -\- V with eigenvector 3~-'^/^(l, 1, 1) and E = 
Ej^ — V/2, which is two fold degenerate with a possible set of orthonormal eigenvectors given by: Vi = 2-i/2(-l, 1, 0) 
and V2 = 6~-^/^(— 1, —1, 2). We note that the remaining two-fold degeneracy is protected by topology and thus persists 
beyond a pcrturbativc description. This degeneracy occurs in the ground state for V > Q (honeycomb lattice) and in 
the first excited state for F < (triangular lattice). 

The remaining two- fold degenaracy of vi,V2 disappears for finite q. The effective two-band Hamiltonian spanned 
by Vi ,V2 is given by: 

(Feff(q))i,- = {vi\H{ci)\vj) (A2) 
iJeff(q) = (i?^ - V/2 + El)l + hvF ||-y + ^qy^ a, + kvp (^q, + ^-q^ a, (A3) 

with vp = 27Th/3am. The eigenenergies are given by i^K+q = E^ ± hvpq with E^ = E^ — V/2 + E'^. We note that 

the effective Hamiltonian can be transformed by a unitary transformation to i?eff = EqI + hvpiqxf^x + qy'^y)^ which 
is identical (except for the value of vf) to the low-energy Hamiltonian obtained in the tight-binding approximation . 




APPENDIX B: TIGHT-BINDING DESCRIPTION FOR SYMMETRIC POTENTIALS 



1. Symmetric honeycomb lattice 



For the honeycomb lattice there are two degenerate minima per unit ceU and for sufRciently deep potential wells 
[i.e. for V <^ 2ma'i ' '^ith the atomic mass and V the confinement amplitudes of Eq. (3)] the basis for the two lowest 
bands can be restricted to the ground state of each minimum. In the symmetric case (equal laser intensities) the 
potential around each minimum is rotationally symmetric and harmonic and the minima form a honeycomb lattice. 
Thus the nearest-neighbor hoppings t are all the same. 

= -t^at(r,)^6(r.+d,)+H.c. = -^(al,,6[) t ) ( J ) ' & = i ^ ^^pH'' ' ^.0 (Bl) 

Here a(r) destroys a particle at the A-site at position r and the dj denotes the vectors connecting an A-site with its 
neighboring B-sites, see Fig. lb. Fourier transformation is given by a{ri) = iV~^/^ exp(ik • ri)ak, where N is the 
number of unit cells. The eigenenergies of the above Hamiltonian are given by: E'k = = ±^[3 + 2cos(fcj,a) + 

4 cos {V3ak^/2)cos{aky/2)]^^^. The energy spectrum is symmetric around E = 0, where the Fermi surface only consists 
of singular (Dirac) points located at the K-points. Close to the K-point one can approximate ^K+q — —vf{(1x — iqy) 
with vp = so that the Hamiltonian close to the K-point can be written as H^s = hvF{<lx<^x + Qy'^y) with 

eigenenergies Eq = zthvpg- 



2. Symmetric triangular lattice 

The triangular lattice has one minimum per unit cell. For a symmetric lattice the potential around the minimum is 
rotationally symmetric and harmonic. The ground state of such a minimum is non-degenerate while the first excited 
state is two fold degenerate and can be represented by Px(x,y) = ipi{x)^po{y) and Py{x,y) = ipo{x)'ipi{y), where x,y 
are measured with respect of the center of the minimum and "00(2;), ipi (x) denote the ground state and the first excited 
state of the one-dimensional harmonic oscillator. 

The band arising from the ground state is described by the following tight-binding Hamiltonian: 

H = -t^clcj =^Ekclck; £;k = -2t[cos(k-ai) + cos(k-a2) + cos(k-(ai-a2))] (B2) 

The hopping between the Px and py orbitals is no longer isotropic^^. As visualized in Fig. 5b, we introduce the 
hopping amplitudes to- and depending on the relative orientation of the orbitals. Non-parallel hoppings (e.g. from 
Px to Py) are excluded by the symmetry of the wavefunctions. 

The tight-binding Hamiltonian describing the two bands arising from the two-fold degenerate first excited state is 
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then given by: 



H 



''k' "k 



Ck 

rfk 



-ffk = /ll 1 + /l2Cra; + h^CTz ; -Ek = /ll ± 



1/2 1/2 
-1 1 



cos(afca;) + 2cos(afcj;/2) cos(afcy\/3/2)j 



-x/Sfr sin(afca:/2) sin(aA:j,-\/3/2) ; hs = U cos{akx/2) cos{akyV^/2) — cos{akx) 



(B3) 



We note that both bands touch at the T point as wcU as at the K-points (see Figs Ic and 2). The energies are given 
by Er = 6tc and Ek = —3tc, while at the M-point the energies are spht with Em = —'^tc ± 2tr- Around the K-points 
k = K + q the energy dispersion to linear order in q is rotationally symmetric, -BK+q = —Stc ± 3^/^atrg/4 + O(q^), 
so that a circular Dirac cone is formed around each K-point. We note that the electronic structure of a triangular 
lattice with degenerate d-orbitals shows similar features"^^. 



4 



* Electronic address: bwunsch@fis.ucm.es 

^ K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 
306, 666 (2004). 

^ A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim (2007), arXiv:0709.1163. 
^ A. K. Geim and K. S. Novoselov, Nature Materials 6, 183 (2007). 
C. Beenakker (2007), arXiv:0710.3848. 

Y. B. Zhang, Y. W. Tan, H. L. Stormer, and P. Kim, Nature 438, 201 (2005). 

K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, 
Nature 438, 197 (2005). 

G. Grynberg, B. Lounis, P. Verkerk, J.-Y. Courtois, and C. Salomon, Phys. Rev. Lett. 70, 2249 (1993). 
S.-L. Zhu, B. Wang, and L.-M. Duan, Phys. Rev. Lett. 98, 150404 (2007). 
C.-Y. Hou, C. Chamon, and C. Mudry, Phys. Rev. Lett. 98, 186809 (2007). 
I. F. Herbut, Phys. Rev. Lett. 99, 206404 (2007). 

W. Zhang, P. Zhang, S. Duan, and X.-G. Zhao (2008), arXiv:0805.2736. 

A. M. Dudarev, R. B. Diener, I. Carusotto, and Q. Niu, Phys. Rev. Lett. 92, 153005 (2004). 
C. Wu and S. Das Sarma (2007), arXiv:0712.4284. 

C. Wu, D. Bergman, L. Balents, and S. Das Sarma, Phys. Rev. Lett. 99, 070401 (2007). 

C. Wu (2008), arXiv:0805.3525. 

L. B. Shao, S.-L. Zlm, L. Slicng, D. Y. Xing, and Z. D. Wang (2008), arXiv:0804.1850. 
F. Haldane, Phys. Rev. Lett. 61, 2015 (1988). 

F. D. M. Haldane and S. Raghu, Phys. Rev. Lett. 100, 013904 (2008). 
S. Raghu and F. D. M. Haldane (2006), arXiv:cond-mat/0602501. 
R. A. Sepkhanov, J. Nilsson, and C. W. J. Beenakker (2008), arXiv:0803.1581. 

J. Mafics, F. Guinea, and A. H. Vozmcdiano, Phys. Rev. B. 75, 155424 (2007). 

J. R. Enshcr, D. S. Jin, M. R. Matthews, C. E. Wicman, and E. A. Cornell, Phys. Rev. Lett. 77, 4984 (1996). 

0. Vafck, Phys. Rev. Lett. 98, 216401 (2007). 

J. Gonzalez, F. Guinea, and V. A. M. Vozmediano, Nucl. Phys. B 424, 595 (1994). 

B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New J. Phys. 8, 318 (2006). 
E. H. Hwang and S. Das Sarma, Phys. Rev. B. 75, 205418 (2007). 
M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, Nature Phys. 2, 620 (2006). 

D. L. Bergman, C. Wu, and L. Balents (2008), arXiv:0803.3628. 
M. V. Berry Proc. R. Soc. Lond. A 392, 45 (1984). 

5. Zhou, G.-H. Gweon, A. Fedorov, P. First, W. de Heer, D.-H. Lee, F. Guinea, A. C. Neto, and A. Lanzara, Nature Materials 

6, 770 (2007). 

31 I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008). 

1. Guillamon, H. Suderow, F. Guinea, and S. Viera, Phys. Rev. B. 77, 134505 (2008). 



32 



